有很多相关的方法都成为MCMC方法。对于网格近似方法而言,我们需要根据给定的点计算先验和似然,并且近似出整个后验分布。MCMC方法要比网格近似更好,这是因为它设计的思想就是将更多的时间用于高概率区域的计算。事实上MCMC方法会根据相对概率访问参数空间中的不同区间。如果区间 $A$ 的概率是区间 $B$ 的两倍,那么我们从区间 $A$ 采样的次数可能是从$B$区间采样次数的两倍,因此即使我们不能从分析的角度得到整个后验,我们也可以采用 $M C M C$ 的方法从中采样, 并且采样数越多, 效果越好。
娶理解 MCMC 方法, 我们先将其拆分成两个 $M C$ 部分, 即蒙特卡洛部分和马尔科夫链部分。
蒙特卡洛这部分可以用随机数的应用来解释。蒙特卡洛方法是一系列应用非常广泛的算法, 其思想是通过随机采样来计算或模拟给定过程。
蒙特卡洛是位于摩纳哥公国的一个非常有名的城市, 蒙特卡洛方法的开发者之一是 Stanislaw Ulam。尽管很多问题都难以求解甚至无法准确用公式表达, 但我们可以通过采样或者模拟来有效地研究。据传, 起因是为了计算单轮游戏中的概率,解决该问题的方法之一是组合分析法; 进行多次单轮游戏, 最后计算其中有多少次是我们感兴趣的。
这听起来似乎是显而易见的, 或者至少是相当合理的, 比如, 你可能已经用重采样的方法来解决统计问題。不过这个实验是早在 70 年前进行的,当时, 第一台计算机才刚开始研发。该方法首次应用于一个核物理问题。如今, 个人计算机都已经足哆强大. 用蒙特卡洛方法可以解决许多有趣的问题,因此, 该方法也被广泛应用到科学、工业以及艺术等各个方面。 可以采用以下下过程估计 $\pi$ :
由此估计得到$\pi$的值
马尔科夫链是一个数学对象, 包含一系列状态以及状态之间的转移概率, 如果每个状态转移到其他状态的概率只与当前状态相关, 那么这个状态链就称为马尔科夫链。
有这样一个马尔科夫链之后, 我们可以任取一个初始点, 然后根据状态转移概率进行随机游走。假设能够找到这样一个马尔科夫链, 其状态转移概率正比于我们想要采样的分布 (如贝叶斯分析中的后验分布), 采样过程就变成了简单地在该状态链上移动的过程。
那么, 如何在不知道后验分布的情况下找到这样的状态链呢? 有一个概念叫做细节平衡条件 (Detailed Balance Condition), 直观上讲, 这个条件是说, 我们需要采用一种可逆的方式移动 (可逆过程是物理学中的一个常见的近似)。也就是说, 从状态 $i$ 转移到状态 $j$ 的概率必须和状态 $j$ 转 移到状态 $i$ 的概率相等。
总的来说就是, 如果我们能够找到满足细节平衡条件的马尔科夫链, 就可以保证从中采样得到的样本来自正确的分布。保证细节平衡的最流行的算法是 Metropolis-Hasting 算法。
为了更形象地理解这个算法, 我们用下面这个例子来类比。
假设我们想知道某个湖的水容量以及这个湖中最深的点, 湖水很浑浊以至于没法通过肉眼来估计深度. 而且这个湖相当大, 网格近似法显然不是个好办法。为了找到一个采样策略, 我们请来了两个好朋友小马和小萌。经过讨论之后想出了如下办法, 我们需要一个船 (当然,也可以是竹筏) 和一个很长的棍子, 这比声呐可便宜多了, 而且我们已经把有限的钱都花在了船上。 (1) 陏机选一个点,然后将船开过去。 (2) 用棍子测量湖的深度 (3) 将移到另一个地点并重新测量。 (4) 按照如下方式比较两个点的结果
- 如果新的地点比旧的地点水位深,就在笔记本上记录新的测量值,并重复过程(2)
- 如果新的地点比旧的地点水位浅,那么我们有两个选择:接受或者拒绝。接受以为记录新的地点并重复(2);拒绝就意味着重新回到上一个点,在此记录上一个点的测量值。
如何决定接受还是拒绝新的测量值呢?这里有一个技巧便是Metropolis-Hasting准则,即接受新的概率等于接受新旧概率值之比。
按照以上过程迭代下去, 我们不仅可以得到整个湖的水容量和最深的点, 而且可以得到整个湖底的近似曲率。在这个类比中, 湖底的曲率其实就是后验分布, 而最深的点就是后验的众数。
迭代的次数越茗, 近似的效果越好。事实上, 理论保证了在这种情形下, 如果我们能采样无数次, 最终能得到完整的后验。幸运地是, 实际上对于很多问题而言, 我们只需要相对较少地采样就足以得到一个相当准确的近似。
现在让我们从更正式的角度来看看该算法。 对于很多分布而言(比如正态分布),我们有相当高效的算法对其采样, 但对于一些其他分布, 情况就变了。 Metropolis-Hastings 算法使得我们能侈从任意分布中以概率 $p(x)$ 得到采样值, 只要我们能筫出某个与 $p(x)$ 成比例的值。 这一点很有用, 因为在类似贝叶斯统计许多问题中,最难的部分是计算归一化因子,也就是贝叶斯定理中的分母。
Metropolis-Hastings 算法的步骤如下。
- 给参数 $x_i$ 赋一个初始值, 通常是随机初始化或者使用某些经验值。
- 从某个简单的分布 $Q\left(x_{i+1} \mid x_i\right)$ 中选一个新的值 $x_{i+1}$ ,如高斯分布或者均匀分存。这一步可以看做是对状态 $x_i$ 的扰动。
- 根据 Metropolis-Hastings 准则计算接受一个新的参数值的概率: $p_a\left(x_{i+1} \mid x_i\right)=\min \left(1, \frac{p\left(x_{i+1}\right) q\left(x_i \mid x_{i+1}\right)}{p\left(x_i\right) q\left(x_{i+1} \mid x_i\right)}\right)$
- 从位于区间 $[0,1]$ 内的均匀分布中随机选一个值,如果第(3)步中得到的概率值比该值大,那么就接受新的值, 否则仍保持原来旳值。
- 回到第 (2) 步重新迭代, 直到我们有足够多的样本, 稍后会解释什么叫足够多。
有几点需要注意。
最后, 我们会得到一连串记录值, 有时候也称采样链或者迹。如果一切都正常进行, 那么这些采样值就是后验的近似。在采样链中出现次数最多的值就是对应后验中最可能的值。该过程的一个优点是: 后验分析很简单,我们把对后验求积分的过程转化成了对采样链所构成的向量求和的过程。
MCMC 方法, 包括 Metropolis-Hastings, 都在理论上保证如果采样次数足够多, 最终会得到后验分布的准确近似。不过, 实际中想要采样足够多次可能需要相当长的时间, 因此, 人们提出了一些 Metropolis-Hastings 算法的替代方案。这些替代方案, 包括 Metropolis-Hastings 算法本身, 最初都是用来解决统计力学中的问题。统计力学是物理学的一个分支, 主要研究原子和分子系统的特性。
汉密尔顿蒙特卡洛方法, 又称混合蒙特卡洛 (Hybrid Monte Carlo, HMC), 是这类改 进方案之一。
简单来说, 汉密尔顿这个词描述的是物理系统的总能量, 而另外个名称中的 “混合” 是指将 Metropolis-Hastings 算法与分子力学 (分子系统中广 泛应用的一种仿真技巧)相结合。HMC 方法本质上和 Metropolis-Hastings 是一 样的, 改进的地方在于: 原来是随机放置小船, 现在有了一个更聪明的办法, 将小船沿着湖底方向放置。为什么这个做法更聪明? 因为这样做避免了 MetropolisHastings 算法的主要问题之一: 探索得太慢而且采样结果自相关 (因为大多数采样结果都被拒绝了)。
那么, 如何才能不必深入其数学细节而理解汉密尔顿蒙特卡洛方法呢? 假设我们还是在湖面上坐着船, 为了决定下一步将要去哪, 我们从当前位置往湖底扔了一个球, 我们假设球面是理想的, 没有摩擦, 因而不会被泥巴和水减速。扔下球之后, 让它滚一小会儿,然后把船划到球所在的位置。现在利用 Metropolis-Hastings 算法中提到的 Metropolis 准则来选择接受或者 拒绝, 重复整个过程一定次数。改进后的过程有更高的概率接受新的位置, 即使它们现在的位置相比前一位置距离较远。
现住挑出我们的思维实验, 回到现实中来。其于汉密尔顿方法䨒要计算函数的梯度。梯度是多维上导数的推导。我们可以利用梯度信息模拟球在曲面上移动的过程。因此, 我们而临一个权衡; HMC计算过程要比 Metropolis-Hasting方法更复杂,但接受效率更高。对于一些复杂的问题HMC要更合适一些。 HMC方法的另一个缺点是:想要得到很好的采样需要指定一些参数。如果手动指定,需要指定一些参数,这需要经验。 PyMC3 中有一个相对较新的不掉向采样算法, 该方法被证实可以有效提升HMC算法效率,同时不必手动调参。
参考资料: